
%% This Function compute the vector of mean inequality constraints theta in the paper
%% Inputs
%set des=1 if Y is descrete
function [mu,mu1,mu2,mu3,mu4]=incons(Y,Z,D,des)
%Estimate the 4 conditional probabilities P(D=d|Z=z) see table 2
P11=sum(Z==1&D==1)/sum(Z==1);
P10=sum(Z==0&D==1)/sum(Z==0);
P01=1-P11;
P00=1-P10;
%Generate  the 4 conditional outcomes  Y|D=d,Z=z
y11=sort(Y(Z==1&D==1));
y01=sort(Y(Z==0&D==1));
y10=sort(Y(Z==1&D==0));
y00=sort(Y(Z==0&D==0)); 
%Estimate E(Y|D=1,Z=0) and E(Y|D=0,Z=1)
my01=mean(y01);
my10=mean(y10); 
%Estimate the proportion of always takers in the subsample D=1, Z=1
q=(P10/P11)*(P10/P11<=1)+(P10/P11>1);
%Estimate the proportion of never takers in the subsample D=0, Z=0
r=(P01/P00)*(P01/P00<=1)+(P01/P00>1);

%Estimate the bounds for the always and never takers if bin=1 we use the correction for descrete outcomes described in Huber and Mellace (2010)
if des==0

y11u=mean(y11(y11>=quantile(y11,1-q)));
y11l=mean(y11(y11<=quantile(y11,q)));
y00u=mean(y00(y00>=quantile(y00,1-r)));
y00l=mean(y00(y00<=quantile(y00,r)));
else 
n11=length(y11);
n00=length(y00); 
y11u=mean(y11(round(n11-q*n11)+(round(n11-q*n11)==0):end));
y11l=mean(y11(1:round(q*n11)+(round(q*n11)==0)));


y00u=mean(y00(round(n00-r*n00)+(round(n00-r*n00)==0):end));
y00l=mean(y00(1:round(r*n00)+(round(r*n00)==0)));    
end

%Estimate the moment inequality constraints
mu1=y11l-my01;
mu2=my01-y11u;
mu3=y00l-my10;
mu4=my10-y00u;
mu=-[mu1;mu2;mu3;mu4];

